library(Seurat)
library(harmony)
library(singleCellTK)
library(data.table)
library(tidyverse)
library(tidyseurat)

# read in data ---------
# use url when possible
fread('CRC-I/data/GSE108989_CRC.TCell.S11138.count.txt.gz') -> mtx

mtx[1:5,1:5]

system.time(
  aggr_mtx <- mtx[, lapply(.SD, sum), .SDcols = -c(1:2), by = symbol]
  )

nona_mtx <- aggr_mtx[!is.na(symbol), ,]

sprs_mtx <- nona_mtx |>
  column_to_rownames('symbol') |>
  as("sparseMatrix")

calc_mtx_sparsity <- function(mtx){
  density <- length(mtx@x) / product(mtx@Dim)
  1 - density
}

# smart-seq2 matrix sparsity is ~84%
calc_mtx_sparsity(sprs_mtx)
# typical 10x sparsity is ~95%

sobj <- sprs_mtx |>
  CreateSeuratObject(min.cells = 3)

sobj

VlnPlot(sobj, c('nCount_RNA','nFeature_RNA'), pt.size = 0)

min(sobj$nFeature_RNA)

glimpse(sobj@meta.data)

sobj$mito.ratio <- sobj |> PercentageFeatureSet('^MT-')

VlnPlot(sobj, c('mito.ratio'))

meta_zhang2018 <- read_delim('CRC-I/data/zhang2018_meta_sample.csv')

sobj <- meta_zhang2018 |>
  mutate(.cell = UniqueCell_ID) |>
  left_join(x = sobj, y = _)

sobj |>
  count(Patient_ID)

ggplot(sobj, aes(Patient_ID, fill = sampleType)) +
  geom_bar()

sobj <- sobj |>
  mutate(tissue = str_extract(sampleType, '^.'),
         tissue = case_match(tissue,
                             'N' ~ 'Normal',
                             'T' ~ 'Tumor',
                             'P' ~ 'Blood'),
         orig.ident = str_c(Patient_ID, sampleType, sep = '_'))

sobj <- sobj |>
  quick_process_seurat()

filter(sobj, time == 'D16') %>%
  ggplot(aes(genotype, fill = cell_type)) +
  scale_fill_viridis_d(option = 'turbo') -> p

p + geom_bar() + NoLegend() -> p1
p + geom_bar(position = 'fill') + ylab('proportion') -> p2

p1 + p2

# plot every cell types in facets
sobj %>%
  filter(time == 'D16') %>%
  group_by(genotype, cell_type) %>%
  tally() -> count_data

count_data %>%
  group_by(genotype) %>%
  summarise(sum(n)) -> sum_genotype

count_data %>%
  rowwise() %>%
  mutate(percent = case_when(
    genotype == 'GG' ~ n/sum_genotype$`sum(n)`[1],
    TRUE ~ n/sum_genotype$`sum(n)`[2]
  )) -> percent_data

ggplot(percent_data, aes(genotype, percent))+
  geom_col(aes(fill = genotype), position = 'dodge') +
  #coord_flip() +
  facet_wrap(vars(cell_type), scales = 'free')

# gene-level analysis ---------
SetIdent(sobj, value = 'genotype') -> sobj

SplitObject(sobj, 'cell_type') -> cell_type_list

enrichR::listEnrichrDbs() -> t

get_path_enrich <- function(fobj){
  Seurat::DEenrichRPlot(fobj,
                        ident.1 = 'RR',
                        ident.2 = 'GG',
                        enrich.database = 'GO_Biological_Process_2021',
                        max.genes = 10000,
                        num.pathway = 24,
                        return.gene.list = TRUE)
}

map(cell_type_list, safely(get_path_enrich)) %>%
  transpose() %>%
  pluck('result') %>%
  compact()-> enrich_result

names(enrich_result) -> enriched_types

transpose(enrich_result) %>%
  pluck('pos') %>%
  set_names(enriched_types) %>%
  compact() -> pos_enrich_result

transpose(enrich_result) %>%
  pluck('neg') %>%
  set_names(enriched_types) %>%
  compact() -> neg_enrich_result

write_rds(enrich_result, '../covid19/results/zhu-liu-enrich-result.rds')

read_rds('../covid19/results/zhu-liu-enrich-result.rds') -> enrich_result
